#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#define ITERATION 25

/*The struct that contains the values of the An, B, Pn and Tn*/
typedef struct{
    mpf_t a;
    mpf_t b;
    mpf_t t;
    mpf_t p;
}Data;

int main(){

    Data D[ITERATION];
    int i;
    int n;

    mpf_set_default_prec(40000016);
    mpf_t PI;
    mpf_t temp;

    for(i=0; i<ITERATION; i++){
        mpf_init_set_d(D[i].a, 0.0);
        mpf_init_set_d(D[i].b, 0.0);
        mpf_init_set_d(D[i].t, 0.0);
        mpf_init_set_d(D[i].p, 0.0);
    }

    mpf_init_set_ui(PI, 0L); // pi =0
    mpf_init_set_ui(temp, 0L); //temp = 0


    mpf_init_set_ui(D[0].a, 1L);

    mpf_sqrt_ui(temp, 2L); // temp = sqrt(2)
    mpf_ui_div(D[0].b, 1L, temp); // b0 = 1/sqrt(2)

    mpf_init_set_d(D[0].t, 0.25); //t0 = 0.25

    mpf_init_set_ui(D[0].p, 1L); // p0 = 1

    /*And let's the magic begin: */
    for(n=0; n<ITERATION-1; n++){
        //Value for A(n+1)
        mpf_add(temp, D[n].a, D[n].b);
        mpf_div_ui(D[n+1].a, temp, 2L);

        //Value for B(n+1)
        mpf_mul(temp, D[n].a, D[n].b);
        mpf_sqrt(D[n+1].b, temp);

        //value for T(n+1)
        mpf_sub(temp, D[n].a, D[n+1].a);
        mpf_pow_ui(temp, temp, 2L);
        mpf_mul(temp, D[n].p, temp);
        mpf_sub(D[n+1].t, D[n].t, temp);

        //value for P(n+1)
        mpf_mul_ui(D[n+1].p, D[n].p, 2L);
        //Clear useless Variables
        mpf_clear(D[n].a);
        mpf_clear(D[n].b);
        mpf_clear(D[n].p);
        mpf_clear(D[n].t);
    }

    //value for PI
    mpf_add(temp, D[n].a, D[n].b);
    mpf_pow_ui(PI, temp, 2L);
    mpf_mul_ui(temp, D[n].t, 4L);
    mpf_div(PI, PI, temp);

    mpf_out_str( stdout, 10, 40000016, PI );

    mpf_clear(PI);
    mpf_clear(temp);

    return 0;
}
